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Drag induced lift in granular media 
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Laboratory experiments and numerical simulation reveal that a submerged intruder dragged hor- 
izontally at constant velocity within a granular medium experiences a lift force whose sign and 
magnitude depend on the intruder shape. Comparing the stress on a flat plate at varied inclination 
angle with the local surface stress on the intruders at regions with the same orientation demonstrates 
that intruder lift forces are well approximated as the sum of contributions from flat-plate elements. 
The plate stress is deduced from the force balance on the flowing media near the plate. 

PACS numbers: 



Objects moved through media experience drag forces 
opposite to the direction of motion and lift forces per- 
pendicular to the direction of motion. The principles 
that govern how object shape and orientation affect these 
forces are well understood in fluids like air and water. 
These principles explain how wings enable flight through 
air and the fins generate thrust in water j]J . 

Lift and drag forces are also generated by movement 
within dry granular media, collections of discrete parti- 
cles that interact through dissipative contact forces. Gen- 
eration and control of these forces while moving within 
granular media is biologically relevant to many desert in- 
habitants that dive into or swim within Q sand. Lift 
forces are also relevant to industrial process such as soil 
tillage 0. 

In granular media, lift and drag forces are not as well 
understood as in fluids; movement probes the complex 
fluid/solid behaviors of dense granular flows While 
progress has been made understanding drag forces in slow 
horizontal and vertical drag and impact Q, there has 
been comparatively little work investigating lift forces. 
Studies have examined lift forces for a partially sub- 
merged vertical rod moving horizontally and for a rotat- 
ing plate 0,0], and the drag force on submerged objects 
with curved surfaces however, the lift forces expe- 
rienced by horizontally translated submerged intruders 
have not been explored. 

Experiment and simulation- Experiment and simula- 
tion were employed to investigate the lift (F z ) and drag 
(F x ) forces on simple shapes during horizontal transla- 
tion in granular media (Fig. 1). In experiment long in- 
truders with different cross-sections were dragged within 
a bed of glass beads with particle diameter (PD) of 
0.32±0.02 cm and density (p) 2.47 g/cm 3 (Fig. 1). Drag- 
ging was performed at a constant speed 10 cm/sec with 
the intruder's vertical mid-point at depth d = 12.5 PD 
and its long axis perpendicular to the direction of mo- 
tion. In experiment, / = 31.3 PD long intruders were 
connected at the midpoint to a force sensor (mounted to 
a linear translation stage) by a stiff stainless steel rod 
of diameter 2 PD. Following the method of 0, forces on 
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FIG. 1. (color online) Lift and drag forces in granular media: 
(a) Schematic of the experimental setup, (b) Lift force as 
a function of depth for the cylinder (•), square rod (■) and 
half cylinder (A). Gray region indicates the depth at which 
forces in (c) were measured, (c) Net force on rods measured 
in experiment («— ) and simulation (■<—). Forces (■<—) on the 
intruder surfaces were measured in simulation and are scaled 
by four for better visibility. 



the connecting rod were determined in separate measure- 
ments and subtracted from F x and F z . The grain bed was 
75 PD wide by 53 PD deep by 75 PD long. The initial 
packing state of the grains was prepared by shaking the 
container moderately in the horizontal direction before 
each run. The volume fraction was determined through 
measurements p, total grain mass (M) , and occupied vol- 
ume (V) to be |f = 0.62 ± 0.01 . 

Simulation employed the soft-sphere Discrete Element 
Method (DEM) [L2| method in which particle-particle 
and particle-intruder contact interactions were deter- 
mined by the normal force F n = fc(5 3 / 2 — GnVnd 1 / 2 and 
the tangential force F s = pF n , where 5 is the "virtual 
overlap" and v n is the normal component of the rel- 
ative velocity. F n comprises a Hertzian contact term 
and a velocity dependent normal dissipation 12(. Con- 



stants k = 2 x 10 b kgs m 
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G n = 15 kgs~ 
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and p = P{ PPy pi} = {0.1, 0.27} represent the hardness, 
viscoelastic constant, and particle-particle and particle- 
intruder friction coefficients respectively 13|- The simu- 
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lated grain bed, bounded by frictionless walls, was 75 PD 
wide by 55 PD deep by 78 PD and consisted of 350,000 
particles in a bi-disperse mixture of equal parts 0.3 and 
0.34 cm diameter spheres. Doubling any dimension of the 
grain bed did not change the forces significantly. The ini- 
tial volume fraction 0.62±0.01 was prepared by randomly 
distributing the particles in the volume corresponding to 
the desired volume fraction and then eliminating parti- 
cle overlap [l4j. Using lower volume fractions reduced 



the magnitude of the forces but gave qualitatively sim- 
ilar results. Forces were nearly independent of speed, 
like other drag studies in the non-inertial regime [l(J, 
(speed 5; 40 cm/s for our study). In simulation the in- 
truder dimensions and intruder speed were matched to 
those in experiment and the forces were averaged over 
the steady-state time interval. 

Shape determines lift- In both experiment and simula- 
tion, F z was sensitive to the cross-section of the intruder. 
As shown in Fig. l(b,c), F z for the half-cylinder was 
downward (opposite the orientation of the curved sur- 
face), while for the two vertically symmetric geometries, 
the full cylinder and square rod, F z was positive with 
magnitude larger for the cylinder than for the square. Ex- 
periments with smaller glass beads (0.3mm) gave similar 
results (see supplementary material). \F Z \ increased with 
intruder depth for all intruders (Fig. lb). The lift mech- 



anism is different from the Brazil nut effect |15| which 
results from agitation of the medium by the container. 

Simulation allowed investigation of the surface stress 
distribution responsible for lift and drag on the intrud- 
ers. For all shapes the surface stress was largest along the 
leading surface (Fig. lc green arrows). Due to the linear 
dependence of granular pressure with depth and the finite 
size of the intruder [f| , local stress increased with depth 
along the flat face of the square (Fig. lc). However, for 
curved intruders (e.g. the cylinder in Fig.lc), the magni- 
tude of local stress was primarily determined by the local 
surface orientation. As the local surface tangent became 
more aligned with the intruder velocity the force magni- 
tude became small, supporting observations that surfaces 
parallel to the direction of motion contribute little to the 
drag force 0] . Since the normal force was larger than the 
frictional force, the direction of the local grain-intruder 
reaction force was nearly opposite to the surface normal 
at all points along the intruder's leading surface. 

The dependence of the forces on the local surface ori- 
entation suggests that decomposing the surfaces into dif- 
ferential area elements and summing the forces on those 
elements may describe the net drag and lift experienced 
by the three shapes studied. A similar decomposition 
was successfully used to calculate net drag and thrust on 
an undulatory sand-swimmer Q in the horizontal plane. 

Plates as differential elements- To determine if the 
forces on curved intruders can be understood by decom- 
posing the shape into flat plate elements we now study 
the stresses on a flat plate with tangent angle a varied 
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FIG. 2. Normal (a), and shear (r), stress on the leading sur- 
face of the cylinder as a function of tangent angle a compared 
to the stresses on a plate with the same a. ■, A, •, and T 

represent a cyHnder , a plate , T cy i inder and r piate respectively. 



between 0° and 180° (e.g. a = 0° is along the direc- 
tion of motion). In simulation, a long (I = 31.3 PD), 
thin (0.1 PD), flat plate of finite width w = 7.94 PD was 
dragged horizontally through the granular medium with 
its center 12.5 PD below the initial surface, and the av- 
erage normal (a) and shear stress (t) on the leading side 
of the plate were measured as a function of a (Fig. 2). 
The stresses were asymmetric about a = 90°, and a in- 
creased rapidly for small a, peaking at a ~ 50°. At 
a « 60°, t changed sign, indicating a reversal in grain 
flow along the surface. We define the effective friction 
ratio on the plate as ji pe ff(a) = t/ct, which is zero at 
a ~ 60° and saturates to the expected magnitude of [i vi 
for a > 135° and a « 0° with opposite signs. Remark- 
ably, along the surface of the cylinder (half-cylinder and 
square rod as well), the stresses approximately matched 
the stresses on the plate oriented at the same angle a 
(Fig. 2) . The stresses on the intruders were corrected by 
considering that the depth of the differential element is 
different from the depth of the center of the cylinder and 
assuming linear dependence of the stress on depth. 
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FIG. 3. (color online) (a) The drag (\f x \, blue) and lift (/ z , 
red) components of the stress on a plate as a function of a in 
granular media (■< and A) as compared to a fluid with with 
Re <C 1 (dashed lines) [T(|. Dash-dot red line is granular 
wedge model (see Force model section). 



The near equality between the stresses on plates and 
local surface regions of intruders with the same orienta- 
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tion implies that F z and F x for a translated rod can be 
approximated by the sum of forces from the correspond- 
ing shape built from infinitesimal plates. Resolving the 
stresses on the plate at orientation a into the lab frame 
(xyz) gives the local lift f z (a) and local drag f x {ce) force 
contributions per unit area on the plate element (Fig. 3). 
F z for different shaped rods results from the integration 
of f z (a) contributions along the intruders leading sur- 
face, corrected by a linear depth term, over the intruders 
infinitesimal surface area dA, e.g. F z = J f z (a)(z/d)dA 
(Fig. 4a). Comparison of this integration over the three 
rod shapes with the measured F z from simulation and 
experiment (Fig. 4b) shows good agreement. The orien- 
tation of the leading surface of the cylindrical intruder 
varies from < a < 180° and the asymmetry of f z (a) 
results in a net positive lift. Positive f z at a = 90° is 
responsible for the small lift on the square rod. For the 
half-cylinder, 90 < a < 180°, and F z is negative. 
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FIG. 4. (color online) (a) F z on intruders calculated by in- 
tegration of / z (a). Red outlines on the shapes indicate the 
leading surface and the hatched area indicates the correspond- 
ing range of a and region of integration, (b) F z calculated by 
integration of f z from simulation (red hatched bar) and model 
(gray hatched bar) compared to direct measurement of F~ in 
experiment (black bar) and simulation (red bar). 

To gain insight into the nature of granular drag and lift 
we compare our results to those from low Reynolds num- 
ber fluids where inertia is negligible and viscous forces 
dominate. f x and f z on a plate at angle a in low Re flu- 
ids are symmetric along the direction of motion a = 90° 
(Fig. 3, dashed curves) while the drag and lift forces 
on the plate in granular media are asymmetric about 
a = 90°. This suggests that a granular model is required 
to understand the origin of lift in granular media. 

Force model- In the quasi-static regime of granular 
flow, the force on an intruder can be determined by 
analyzin g th e force balance on the moving volume of 
grains jlHH • With the plate acting as one flow boundary 
we can determine the normal and tangential components 
of stress (tr and r) on the plate surface from these force 
balance equations. In practice this requires approximat- 
ing the boundaries of the moving media as planes and 
computing the forces acting on them. 

Examination of the motion of grains in the vertical 
plane (xz) reveals that the particles move upwards in 



FIG. 5. (color online) Flow of grains and force-balance model: 
(a) Flow held in the vertical (xz) plane for three plate (solid 
black line) orientations. Gray boundary indicates regions 
with upward flow [18]. (b) Forces on a wedge for a < a vp 
[regime (VP)] and on the band for a > a vp [regime (P)]. (c) 
The weight of the upward-flow region as function of a cal- 
culated from simulation (■) and fit (black) to W = csin(a), 
where c = 5.7 N. (d) The average flow velocity angle tp in 
the upward-flow region vs. a (•). (e) Normal component of 
the stress on the plate a calculated from the model (black) 
and measured from simulation (a). The black dashed line 
a — a vp = 97° indicates the boundary between the two 
regimes (VP) and (P). 



front of the plate and flow along a lower boundary (see 
Fig. 4a), a slip plane. Finite yield stress in granular me- 
dia results in flowing regions bounded by slip planes with 
upward flow direction due to increasing yield stress with 
depth [loj • The upper region of the flow is confined by a 
boundary starting from the top edge of the plate and ap- 
proximately parallel to the lower boundary. The weight 
W of this up- flow region A 17J is well fit by W = csin a, 
and is thus proportional to the projected plate length 
normal to the direction of motion (see Fig. 5c). The di- 
rection of average velocity of the particles in the band tp 
varies little (see Fig. 5d) and we approximate the angle 
of the flow boundaries as ip = 44°. 



The plate defines one of the boundaries of A for large 
a (e.g. a = 150° in Fig. 5a) but for small a (e.g. a = 50° 
in Fig. 5a), the vertical velocity of the particles adjacent 
to the plate is negative and thus the upward flow bound- 
ary is described by a virtual plane intersecting the top 
of the plate and extending downwards at angle a vp . We 
therefore consider two regimes of flow: for a > a vp , we 
approximate A as bounded by the plate and two paral- 
lel surfaces with angle -0 [regime (P) in Fig. 5b]. For 
a < ctvp, the upwards flow region A is bounded by the 
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virtual plane and two parallel surfaces with angle ip. The 
region adjacent to the plate is bounded by the plate, the 
virtual plane, and an approximately horizontal bottom 
surface [regime (VP) in Fig. 5b]. For regime (P), the 
forces on the flowing band are the forces from the top 
and bottom boundaries and its weight. On the plate 
surface, the friction coefficient /x pe //(a) is used. Within 
the media, dynamic friction is assumed and the friction 
coefficient is tan 7, determined by the angle of repose, 
7 = 13±1° , which is measured in separate simulations 
by tilting the initially horizontal container and recording 
the post avalanche surface orientation. Simulation indi- 
cates that the stress at the bottom surface dominates; 
therefore we neglect the force on the top of the flow band 
(F t in Fig. 5b). The normal stress on the plate, a, can 
be solved from the force balance on the band to obtain 

= T^ J°(S-t- 7 7) ' where /3(a) = tan ~ 1 ^ff- 
For regime (VP), the normal stress on the plate is 

calculated by considering the stress on the wedge ad- 
jacent to the plate. The stress on the virtual plane can 
be solved by the above equation with the corresponding 
weight of the band and a fixed angle a = a vp . Solv- 
ing the force balance equation for this triangular wedge 
we obtain a (a) = K J^lsfc) ESi^f^l whe re 

V I wl Bm(a vp — f)'— i>— 7) sin(a-/3-7) ' 

P' = P(a vp ). The parameter a vp = 97°, determined from 
a best fit of a(a) from simulation, is within the expected 
range from flow field observations (Fig. 5a). 

Comparison of a calculated from the wedge model and 
a measured directly in simulation demonstrates that the 
model captures the asymmetric shape of cr (Fig. 5e) . In- 
tegration of fx calculated from the model over the non- 
planer intruder surfaces yields net lift forces in agreement 
with those from other methods (Fig. 4b). The decrease 
in cr above a = 90° results from the decrease in W with 
increasing a. For a < 90°, although W increases with 
a, the stress on A is transmitted by the wedge which in- 
duces extra resistance on the bottom plane; therefore a 
peaks at a smaller than 90° . The model assumes the up- 
flow area A (W) increases with depth which explains the 
monotonic increase of \F Z \ with depth. The discrepancy 
between a in model and simulation may be due to the 
simplified description of the shape of the boundaries of 
the flowing media, approximation of ip and a vp as con- 
stants and W as a simple function, and neglecting F t . 

Decomposition of the force on the intruder into forces 
on differential elements assumes that the flowing region 
corresponding to each element is not disturbed by other 
elements. This may explain the difference we observe be- 
tween the stress on the plate and the local stress along 
the intruder for small a (corresponding to elements on 
the bottom), where the upper flow boundary for these 
elements is obstructed by higher regions of the cylinder. 
For objects with concave leading surfaces, the use of dif- 
ferential surface elements may require consideration of 
material jamming in the concave region. 

We have shown that the magnitude and sign of the 



drag induced lift force in granular media depends on the 
shape and depth of the intruder. Drag induced lift on 
non-planar intruders can be computed as the summation 
of lift forces from planar elements which each experience 
a lift force resulting from the pushing of material up a 
granular slip plane. The increase of yield stress with 
depth in granular media causes the asymmetric flow and 
enhances the lift force on a plate facing downward. Our 
understanding of these forces can elucidate the effects 
of head and body shapes [HI of sand burrowing organ- 
isms, and could aid design of control surfaces to allow 
robots [l^l to maneuver in granular environments. 
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